Effect of acetylcholine deficiency on neural oscillation in a brainstem-thalamus-cortex neurocomputational model related with Alzheimer’s disease

Previous works imply that involving brainstem in neuropathological studies of Alzheimer’s disease (AD) is of clinically significant. This work constructs a comprehensive neural mass model for cholinergic neuropathogenesis that involves brainstem, thalamus and cortex, wherein how acetylcholine deficiency in AD affects neural oscillation of the model output is systematically explored from the perspective of neurocomputation. By decreasing synapse connectivity parameters in direct cholinergic pathway from brainstem to thalamus or in indirect glutamatergic synapse pathway from cortex to brainstem to mimic the pathological condition of reduced acetylcholine release in patients with AD, the property of neural oscillation in this model is numerically investigated by means of power spectrum in frequency domain and amplitude distribution in time domain. Simulated results demonstrate that decreasing synapse connectivity whether in the direct cholinergic pathway or in the indirect glutamatergic synapse pathway can alter the neural oscillation significantly in three aspects: it induces an obvious decrease of dominant frequency; it leads to a degraded rhythmic activity in the alpha frequency band as well as an enhanced rhythmic activity in the theta frequency band; it results in reduced oscillation amplitude of the model output. These results are agreement with the characteristic of electrophysiological EEG measurement recorded in AD, especially support the hypothesis that cholinergic deficiency is a promising pathophysiological origin of EEG slowing in AD. Our analysis indicates that targeting the cholinergic system may have potential prospects in early diagnosis and treatment of AD.

www.nature.com/scientificreports/ M2 receptor subtype [9][10][11]43 . Meanwhile, the SC population in the brainstem sends excitatory outputs to the TRC and IN populations via glutamatergic neurotransmitter 40,41 . The synaptic interconnections between the thalamus and the cortex modules conform to the classical results as presented in the literatures 26,36,37 , i.e., all the populations of TRC, IN and TRN in the thalamus region receive excitatory inputs from the PY population, on the same time, the TRC population sends excitatory outputs to all the populations of PY, eIN, fIN and sIN in the cortical region. In addition to the ascending synapses from the brainstem to the thalamus and cortex, there is important descending synapse pathway from the cortex to the brainstem in the central nervous system 44,45 . Here, the PY population in the cortex module sends excitatory outputs to the SC and PNN populations in the brainstem module. Other extrinsic sources to the BTC model come from the retinal and nearby cortical formations, i.e., the retinal population sends excitatory outputs to the IN, TRC and SC populations, and the PY population receives excitatory afferent from nearby cortical regions 26,28,29,36,37 . For convenience, neuron populations in this model are denoted by different lowercase letters s , d , i , t , n , h , p , f , l , r and c , which represent populations SC, PNN, IN, TRC, TRN, eIN, PY, fIN, sIN, Ret, Cor, respectively. According to the technique of neural mass model, the kinetic equations governing each population are obtained by performing two mathematical operations. One is that the average membrane potential V a in each population a (here a = s , d , i , t , n , h , p , f , l , r , c ) receiving from all afferent neuronal populations is changed into an average density of spikes ϕ a , which is usually simulated by the following sigmoidal function: where C abe represents excitatory synaptic connection parameter between postsynaptic neuron population a and presynaptic neuron population b , and C adi represents inhibitory synaptic connection parameter between postsynaptic neuron population a and presynaptic neuron population d . x b and x d are postsynaptic potential (PSP) respectively generated by population b and d . e 0 determines the maximum firing rate, V 0 is the firing threshold, σ controls the steepness of this sigmoid function.Note that the average density of spikes for neuron population Ret is simulated by Gaussian white noise P 1 (t) with mean µ 1 and variance ε 1 . The average density of spikes for neuron population Cor is simulated by Gaussian white noise P 2 (t) with mean µ 2 and variance ε 2 .The other mathematical operation is that a second linear transform function of pulse response converts the presynaptic average spikes density ϕ a into the postsynaptic membrane potential x a , which is as follows: where A a is the synaptic strength determining the maximum amplitude of PSP, and τ a represents the time constant of PSP. Please refer to the detailed mathematical equations governing all the 11 neuron populations in this BTC model in the Supplementary Material section.
The connectivity parameters of afferent to the brainstem and thalamus modules are according to the previous physiological data [46][47][48][49][50] . In detail, based on the electron microscopic observances on the connectivity patterns of two main cell types in the lateral geniculate nucleus of the cat, Erisir et al. 46 have revealed that distribution of synapse is different between relay cells and interneurons present in this nucleus. The relative distribution of terminal types contacting relay cells is 14.6 ± 0.6% retinal (RLP) terminals, 29.6 ± 1.6% inhibitory (F) terminals, and 55.8  www.nature.com/scientificreports/ ± 1.7% cortex and brainstem terminals. Whereas, the analysis of the total synaptic inputs onto the interneurons indicates that interneurons receive 37.8 ± 1.0% of all synapses from retinal terminals, 26.8 ± 1.5% from inhibitory terminals, and 35.4 ± 1.8% from cortex and brainstem terminals. Moreover, in geniculate A-laminae the relay cells and the interneurons receive 85.8 ± 0.6% and 14.2 ± 0.6% of all synaptic terminals, respectively, i.e., the synaptic terminals contacting the relay cells is about 6 times of that contacting the interneurons. Erisir et al. 47 have further found that the relative number of cortical inputs and brainstem inputs to the lateral geniculate nucleus is of the same order, each of whose terminals constitutes roughly one-half. As for the TRN afferent axons, Jones have reported that almost 70% synaptic inputs to the reticular nucleus in the somatosensory sector of the rat is attributed to corticothalamic terminals, 20-25% is attributed to thalamocortical collateral synapses, 15-20% is attributed to GABA ergic synapses and a quite small number of inputs is from the brainstem 48 . Another study has reported that the synaptic proportions from corticothalamic, thalamocortical and GABA ergic synapse are respectively 60-65%, 20% and 15% 49 50 . On the basis of the above studies 46-50 the connectivity parameters relative with the brainstem and thalamus modules are determined. On the same time, the connectivity parameters afferent to the cortex module are sourced from the works 26,28,33,34,51 . The detailed values for all the connectivity parameters are described in Table 1 in the Supplementary Material section. In addition, the synaptic strength A a , time constant τ a , as well as other basic parameters ( e 0 , V 0 , σ , µ 1 , ε 1 , µ 2 , ε 2 ) in the BTC model are referenced from previous works 26,33,37 , which are outlined in Table 2 in the Supplementary Material section.

Main result
As described in the Introduction, previous reports about reduced ACh release, loss of cholinergic neurons and decreased ChAT activity in AD patients have indicated that acetylcholine release is abnormal in AD patients [5][6][7][8] .
In brainstem there are a large number of cholinergic neurons projecting to thalamus, moreover, the brainstem nuclei such as superior colliculus and cholinergic nucleus are susceptible to AD [9][10][11] . This implies that AD is associated with degeneration of the direct cholinergic pathway from brainstem to thalamus. In addition, the neuroimaging technique has shown that loss of functional/structural connectivity in cortical areas have been related to cortical atrophy in AD from graph theoretical studies 52,53 . This property of cortical atrophy may influence the descending synapse projection from cortex to brainstem 44,45 , which in turn could weaken the cholinergic pathway from brainstem to thalamus indirectly. Thus, this study simulates the pathological condition of reduced acetylcholine release in AD by decreasing the synapse connectivity parameters in the direct cholinergic pathway from brainstem to thalamus (i.e., C tde from PNN to TRC, C idi from PNN to IN, C ndi from PNN to TRN) and the indirect glutamatergic synapse pathway from cortex to brainstem(i.e.,C dpe from PY to PNN). How the abnormal acetylcholine release affects neural oscillation in this BTC model is systematically explored by means of power spectrum and amplitude distribution. In the following, the second-order differential equations governing the BTC model are numerically solved by Euler method in an environment of MATLAB 2019a. The total simulation time is 120 s with a time resolution of 1/256 s. The model output corresponds to the summated postsynaptic potential V t in the TRC population. For each output vector, an epoch is extracted in the interval of [20,120] s to ensure that all the transients are discarded. The epoch of output is then bandpass filtered by a Butterworth filter of order 10 with a lower and upper cut-off frequencies of 0.5 and 50 Hz, respectively. The detailed analysis of power spectrum and oscillation amplitude are further carried out based on the filtered output.
Decreasing the excitatory cholinergic connectivity from PNN to TRC . The level of cholinergic neurotransmitter from the PNN population in the brainstem module to the TRC population in the thalamus module is determined by the excitatory synapse connectivity parameter C tde . In this section, through decreasing the strength of C tde to mimic the degeneration of acetylcholine release in patients with AD, neural oscillation of this BTC model is firstly characterized by means of power spectrum analysis in frequency-domain. Based on the filtered output of V t , the power spectrum density of the BTC model is computed using the Welch technique with a Hamming window. Quantitative analysis of power spectrum is performed by characterizing dominant frequency and measuring relative power in specific frequency bands. Dominant frequency is a frequency at which the power spectrum density reaches its peak. The relative power of alpha band (8)(9)(10)(11)(12)(13) or theta band (4-8 Hz) is calculated by averaging the relative power spectrum density within the corresponding frequency band.
The dominant frequency of the model output is illustrated in Fig. 2a when the connectivity parameter C tde is varied in the rage of[60, 100]. On the whole, the dominant frequency shows a downward trend from 9.25 to 6.25 upon decreasing C tde from 100 to 60: it is initially within the alpha band, then the dominant frequency decreases steadily until C tde is decreased to a critical value of about 80, at which the dominant frequency transits from the alpha band to the theta band, afterwards it decreases continually within the theta band as C tde is further decreased. That is to say, the smaller the excitatory cholinergic synapse connectivity parameter C tde , the lower the dominant frequency. This phenomenon is further vividly confirmed by some individual plots of power spectrum density. Figure 2b depicts the detailed power spectrum density when the connectivity parameter C tde is successively chosen as 100, 82, 70 and 60. Clearly, peaks of the power spectrum density in the upper two panels of Fig. 2b are concentrated at the alpha frequency band of 9.25 and 8.5 respectively when C tde is more than 80, whereas peaks in the lower two panels of Fig. 2b are centered at the theta frequency band of 7 and 6 respectively when C tde is less than 80. These interesting results indicate that acetylcholine deficiency, a biomarker for AD, can induce an obvious decrease of dominant frequency, which is consistent with the electrophysiological experiment results that the dominant/peak frequency of EEG is significantly lower in early stage of AD than that in control subjects 15 www.nature.com/scientificreports/ Furthermore, the power spectrum analysis of relative power in specific frequency bands is carried out. The relative power within the alpha band and within the theta band is respectively illustrated in Fig. 2c when the excitatory synapse connectivity parameter C tde is in the range of[60, 100]. One can observe that upon decreasing C tde from 100, the relative alpha band power decreases slightly until C tde arrives at the critical value of about 80, then it falls sharply till C tde ≈ 73 , after which the relative power in alpha band does not decrease anymore and basically tends to be stable. Interestingly, the evolution of the relative theta band power is opposite to that of the relative alpha band power. In detail, upon decreasing C tde from 100 to 60, the relative power in theta band is just a little increased till C tde reaches the critical value of about 80, then it grows swiftly until C tde ≈ 73 , after that the relative theta band power fluctuates slightly with the further decrease of C tde . This phenomenon accords with the traditional experiment EEG measurements that there is a decrease in the alpha band activity and an increase in theta band activity in patients with AD 17-21 , in particular, it supports the hypothesis that cholinergic deficiency is a promising pathophysiological origin of the EEG slowing in AD 20,22 .
Secondly, neural oscillation of this BTC model is investigated by means of oscillation amplitude in timedomain. When decreasing the excitatory synapse connectivity parameter C tde from 100 to 60, the model output, i.e., the summated postsynaptic potential V t in the TRC population, is illustrated in Fig. 3a for each value of C tde . Obviously, as a whole the oscillation amplitude of V t is getting smaller upon decreasing the synapse connectivity parameter C tde . To visualize this result, the detailed postsynaptic potential V t for some typical synaptic strengths such as C tde = 100, 82, 70 and 60 is also depicted in Fig. 3b. one can qualitatively observe that the oscillation amplitude of V t is gradually reduced with the decrease of connectivity parameter. Especially, for the two former ones the BTC model produces rhythmic oscillation of potential analogous to alpha-rhythm (with the dominant frequency within alpha band as shown in Fig. 2b), i.e., they exhibit waxes and wanes in amplitude just like a form of spindle 32 . There is no obvious qualitative change in dynamics from the perspective of time series. As the synapse connectivity C tde is less than the critical value of 80, the output of BTC model in the latter two panels presents rhythmic oscillation of the theta frequency band (with the dominant frequency within the theta band as shown in Fig. 2b) with relative small amplitude compared to the former two cases. Due to the influence of the noise, statistical and uncertainty analysis is performed in order to avoid the randomness. Furthermore, as displayed in the Fig. 3c the statistical property of oscillatory potential is characterized by amplitude distribution based on 500 realizations of postsynaptic potential. Here the histogram of oscillation amplitude is fitted by normal density function with estimated parameters of mean and variance. As can be seen in the Fig. 3c, the estimated mean of oscillation amplitude is 11.04, 8.68, 5.87 and 3.43 when the connectivity parameter C tde is successively taken as 100, 82, 70 and 60. This result further quantitatively confirms that the reduced oscillation amplitude of the model output results from the decreased excitatory synapse connectivity parameter C tde , i.e., a hallmark of deficit cholinergic synapse pathway from the PNN population to the TRC population. www.nature.com/scientificreports/ Decreasing the inhibitory cholinergic connectivity from PNN to IN and TRN. In addition to the above excitatory cholinergic pathway from PNN to TRC via M1 and M3 muscarinic receptor subtypes, there are also two inhibitory cholinergic pathways from PNN to IN and TRN via M2 muscarinic receptor subtype, whose synaptic connection are determined by C idi and C ndi , respectively. In this section, the physiological feature of acetylcholine deficiency in AD patients is reflected by reducing C idi or C ndi . Under such environment, the neural oscillation in the BTC model is investigated by means of power spectrum in frequency domain and amplitude distribution in time domain. Firstly, neural oscillation of this BTC model is discussed on the basis of power spectrum analysis. The variation of dominant frequency of the model output V t with the decrease of inhibitory cholinergic synaptic strength is delineated in Fig. 4a and b. One can see that whether decreasing the synapse connectivity parameter C idi (from PNN to IN) or C ndi (from PNN to TRN), they can always lead to a decrease of dominant frequency. To be specific, when C idi (or C ndi ) is decreased from 18 (or 10) to a critical value of about 5 (or 2), the dominant frequency gradually drops within the alpha band, then it steps into the theta band and continually drops within theta band. To illustrate the above result more detail, some power spectrum density curves are exhibited for different synaptic connectivity parameter values. In Fig. 4c, C idi is taken as 10, 6, 3 and 1 from left to right. Clearly, when C idi is greater than 5, the frequency corresponding to the maximum of power spectrum density is respectively at 9.25 and 8.25, whereas when C idi is less than 5, the frequency corresponding to the maximum of power spectrum density is respectively at 7 and 6. That is to say, the dominant frequency for the left two figures and the right two figures is within the alpha band and theta band respectively. In Fig. 4d, C ndi is taken as 7, 3, 1.5 and 0.5 from left to right. Obviously, when C ndi is larger than 2, the dominant frequency in the left two panels is at the alpha frequency band of 9.25 and 8.5 respectively, while the dominant frequency in the right two panels is at the theta frequency band of 7.5 and 7.25 respectively when C ndi is smaller than 2. In addition, the relative power of the model output is calculated with the change of inhibitory cholinergic synapse connectivity C idi and C ndi in Fig. 4e and f. From the two panels, we can observe that upon decreasing synapse connectivity the relative theta band power (denoted by red curves with circles) increases mildly until C idi (or C ndi ) reaches the critical value of about 5 (or 2), then it increases sharply with the further decrease of synapse connectivity. Interestingly, the relative alpha band power (denoted by blue curves with diamonds) decreases slowly till C idi (or C ndi ) reaches the critical value of about 5 (or 2), after that it decreases swiftly with a further reduction of synapse connectivity.
In a word, the reduction of inhibitory cholinergic synaptic strength from the brainstem to the thalamus can change the rhythm of neural oscillation in the BTC model significantly. It can induce not only a diminished dominant frequency but also a slowing of rhythmic content, i.e., a decreased alpha band activity and an increased theta band activity. These simulated phenomena resulting from acetylcholine deficiency agree with the EEG characteristics observed in the electrophysiological experiments of AD patients, i.e., the major effect of AD on EEG is slowing of EEG 18-22 along with the degraded peak frequency [15][16][17] .
Secondly, the neural oscillation of the BTC model is explored in time domain by analyzing oscillation amplitude of model output. Figure 5a and b depict the variation of the summated postsynaptic potential V t during [20,120]s as the inhibitory cholinergic synapse connectivity is diminished. These two figures indicate that the oscillation amplitude is degraded as a whole upon diminishing the connectivity parameter C idi or C ndi , though this downward trend is not very smooth. In order to confirm this phenomenon, the detailed amplitude distribution at some synaptic strengths is further displayed in Fig. 5c and d according to the statistical property of 500 realizations of postsynaptic potential V t . From the fitted normal density function of oscillation amplitude, the Fig. 5c reveals that the estimated mean of oscillation amplitude is successively 9.58, 8.66, 6.92 and 5.41 when C idi = 10, 6, 3 and 1. On the same time, the Fig. 5d reveals that the estimated mean of oscillation amplitude decreases from 9.83 to 8.82, then to 7.76, and finally to 6.89 when C ndi decreases from 7 to 3, then to 1.5, and finally to 0.5. These results imply that decrease of inhibitory cholinergic synapse connectivity from brainstem to thalamus resulting from acetylcholine deficit can indeed bring about a decrease in oscillation amplitude of the model output.
Decreasing the excitatory glutamatergic connectivity from PY to PNN. As illustrated by the schematic of BTC model in Fig. 1, the release of acetylcholine from brainstem to thalamus is indirectly modulated by glutamatergic synapse pathway from cortex to brainstem. The synapse loss caused by cortical atrophy in AD may destroy the descending synapse projection from cortex to brainstem 44,45 . In this section, decreasing the glutamatergic synapse connectivity C dpe from PY in cortex to PNN in brainstem to mimic indirect acetylcholine deficiency, neural oscillation in this BTC model is discussed using power spectrum analysis and oscillation amplitude.
Firstly, the property of neural oscillation is characterized by the dominant frequency and the relative power in frequency domain. The evolution of dominant frequency of the model output V t during the decrease of C dpe is depicted in Fig. 6a. Clearly, as C dpe is diminished from 78 to 38, the dominant frequency initially falls within the alpha frequency band until C dpe arrives at a critical value of 45, and then the dominant frequency continues to fall within the theta frequency band with a further decrease of C dpe . Individual curves of power spectrum density at some synaptic connectivity are further given in Fig. 6b to confirm the above phenomenon. As shown in the upper two panels ( C dpe = 65, 48), the peak value of power spectrum density is obtained at about 9.25 and 8.5, i.e., the dominant frequency is located within the alpha frequency band when the synaptic connectivity is large than the critical value of 45. Nevertheless, as displayed in the lower two panels ( C dpe = 41, 38), the dominant frequency is respectively 7 and 6.5, which indicates the peak of power spectrum density is located within theta frequency band when the synaptic connectivity is smaller than the critical value of 45. Furthermore, the variation of relative power along with the synaptic connectivity C dpe is displayed in Fig. 6c. Upon decreasing C dpe from 78, the relative alpha band power in the left panel fluctuates slightly until C dpe reaches a certain value of 45, then it falls sharply till C dpe ≈ 38 . Interestingly, an opposite phenomenon occurs for the evolution of the www.nature.com/scientificreports/ www.nature.com/scientificreports/ relative theta band power in the right panel, i.e., the relative power within theta band rises slightly until C dpe arrives at the certain value of 45, then it booms quickly till C dpe ≈ 38 . As we expected, the rhythmic property of neural oscillation in this BTC model resulting from glutamatergic synapse deficit is similar with that what happened in the case of direct acetylcholine deficiency, i.e., the diminished dominant frequency, the degraded alpha band activity together with enhanced theta band activity, which is consistent with the EEG characteristic obtained from clinical trials of AD patients 15,17,20,22 . Secondly, neural oscillation of model output in the BTC model is studied by means of oscillation amplitude in time domain. Figure 7a depicts the summated postsynaptic potential V t of model output during [20,120] s of simulations for every connectivity parameter C dpe ranged in [38,78]. Intuitively, as C dpe is decreased from 78 to 38 the oscillation amplitude of V t slowly descends as a whole except for few connectivity parameter at which the oscillation amplitude is suddenly enlarged. In addition, this result is visualized by some individual postsynaptic potential V t at different synaptic strength such as C dpe = 38, 41, 48, and 65. Combined with the dominant frequency shown in Fig. 6a, the left two panels in Fig. 7b reveal that the model output has an visible alpha rhythmic content with the amplitude waxing and waning when the synapse connectivity is greater than the critical value of 45. Whereas, in the right two panels of Fig. 7b, there appears an obvious theta rhythmic content in the postsynaptic potential with relatively small oscillation amplitude when the synapse connectivity is less than the critical value of 45. Furthermore, the corresponding amplitude distribution of 500 realizations of postsynaptic potential is illustrated in Fig. 7c. One can observe that the estimated mean of oscillation amplitude is in turn 9.63, 8.96, 6.57 and 4.83 on the basis of the fitted normal density function of oscillation amplitude. This quantitative result once again confirms that the decreased oscillation amplitude of model output is due to the damaged glutamatergic synapse connectivity C dpe from PY to PNN, i.e., an indirect acetylcholine synapse pathway from cortex to brainstem.

Conclusion and discussion
This work establishes a comprehensive neural mass model in the interactive brain structures of brainstem, thalamus and cortex. It takes destroying cholinergic synapse pathway as a surrogate for the abnormal cholinergic system in patients with AD. By monitoring synapse connectivity strength in direct cholinergic pathway from brainstem to thalamus and indirect glutamatergic synapse pathway from cortex to brainstem, we are able to simulate some property of neural oscillation in this model under an environment of acetylcholine deficiency from the standpoint of neurocomputation. By analyzing power spectrum of the model output in frequency domain, the results reveal that upon diminishing synapse connectivity strength in a certain range the dominant frequency initially decreases within alpha frequency band and then steps into theta frequency band, which is consistent with the electrophysiological experiment results that the dominant/peak frequency of EEG is significantly lower in early stage of AD than that in control subjects [15][16][17] . Meanwhile, the neural oscillation presents a slowing rhythmic content with a degraded alpha band relative power and an enhanced theta band relative power, which accords with the traditional experiment EEG measurements that there is a decrease in the alpha band activity and an increase in the theta band activity in patients with AD 17-21 . In particular, it supports the hypothesis that cholinergic deficiency is a promising pathophysiological origin of the EEG slowing in AD 20,22 . What's more, amplitude www.nature.com/scientificreports/ distribution of the neural oscillation in time domain suggests that the oscillation amplitude is overall diminished with the reduction of synapse connectivity. We expect these findings could have important implications on better understanding cholinergic pathogenesis and expounding potential feature for AD. At last, we point out that this work, by reducing the synaptic connection parameters in cholinergic direct pathway from brainstem to thalamus or indirect glutamatergic synapse pathway from cortex to brainstem, mainly focus on the effect of acetylcholine deficiency on brain neuron activity in patients with AD. As known that there could be neuron loss and brain atrophy in cortical region in AD patients 44,45 , which would induce unfavorable synapse information processing. Thus, reducing the synaptic connection parameters relative to cortical region can also simulate the pathological state of AD. Through further simulation, we find that when synapse connectivity parameter decreases from PY to TRC, the change of the dominant frequency and the oscillation amplitude is similar with what reported in this work, i. e. there is an obvious decrease of dominant frequency, a degraded rhythmic activity in the alpha frequency band as well as an enhanced rhythmic activity in the theta frequency band, and the reduced oscillation amplitude of the model output. For the limited space, we do not list them in detail.

Data availability
All data generated or analyzed during this study are included in this published article.